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Abstract. - We study gravitational clustering of mass points in three dimensions with random 
initial positions and periodic boundary conditions (no expansion) by numerical simulations. 
Correlation properties are well defined in the system and a sort of thermodynamic limit can be 
defined for the transient regime of clustering. Structure formation proceeds along two paths: 
(i) fluid-like evolution of density perturbations at large scales and (ii) shift of the granular (non 
fluid) properties from small to large scales. The latter mechanism finally dominates at all scales 
and it is responsible for the self-similar characteristics of the clustering. 



One of the fundamental challenges of modern cosmology is the understanding of the for- 
mation of the structures in the Universe. These structures consist of clusters of galaxies and 
show complex properties extended to very large scales ffl. Usually the simulations and the 
models aimed at the understanding of these structures are based on three essential elements: 
(i) the dynamics under the effect of the gravitational forces; (ii) some particular type of initial 
conditions; (iii) a model for the Hubble expansion |^,|^]. In addition simulations are usually 
run up to a time which is supposed to represent the present state. 

Here we would like to take inspiration from these studies and formulate the problem 
of clustering by gravity in the perspective of the statistical physics of dynamical systems. 
So we will single out the role of each individual effect at the expense of a loss in realism. 
We try therefore to identify simple fundamental mechanisms which can be studied in great 
detail and followed up to their asymptotic state. As a first problem we consider the simple, 
basic question: how does a random distribution of point masses evolve under gravity? The 
comparison with the expanding case may then allow us to identify the specific role of this 
effect. Simulations similar to ours, but in a cosmological context, were performed long ago, 
e.g., by Itoh et al. Q. However, we are going to see that the general problematic we consider 
S_j ' and the final interpretation will be substantially different. 

The main results are: (i) the existence of a well defined thermodynamic limit for the 
transient regime; (ii) the nature of the clustering process arising from the shift of the granular 
(non fluid-like) characteristics from small to large scales, (iii) the evolution of correlations 
shows self similar characteristics. 
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Self-gravitating systems have been studied in various perspectives since a long time. How- 
ever, due to the peculiarities of the 1/r potential, most of the concepts and methods of 
standard statistical physics are problematic Indeed, the usual approach of statistical 

mechanics cannot be applied to a system of many point particles interacting by the Newtonian 
potential, because of (i) the long range nature of the 1 /r potential and of (ii) the divergence at 
the origin ||. In particular, energy is not extensive, therefore the thermodynamic properties 
cannot be uniquely defined as the microcanonical and canonical ensembles are not equivalent. 

In order to simulate an infinite system, N particles at rest are placed at random into 
a cube of side L with periodic boundary conditions. Every particle in the simulation box 
interacts with all the other particles and with the periodic "replicas" of the whole system. 
The Hamiltonian for a system of N identical particles of mass m is: 

* (ij) 

where the second summation runs over all pairs. The interaction potential which accounts for 
both the pair contribution and that of the replicas, is: 

$ (r) = > h / dr 

U ^ [ \r + Ln\ + J V |r' + Ln| a . 

where G is the gravitational constant, V is the box domain and n runs over all three di- 
mensional vectors with integer components. In the above integral the first term is the usual 
Newtonian potential due to a particle with mass m placed in the replica specified by n, the 
second term is introduced to remove the long range divergence of the potential. It can be 
seen as the potential generated by a repulsive mass —m smeared uniformly over the volume 
of the same replica. The summation over the replicas is given by the summation of n. The 
system is individuated by n = 0. Such definition for 4> L gives the expected expression for 
the force due to an infinite periodic system. Periodic boundary conditions have two major 
advantages with respect to free boundaries, when simulating an infinite system with a finite 
N-body representation, namely: (i) they avoid having a preferential point in the simulation 
and (ii) the total number of particles inside the cube of size L is constant, so particles cannot 
"evaporate" , as would happen in an open system Jl0[| . Furthermore they ensure that both 
terms in Eq. (Q) grow linearly with N for constant density. In order to simulate the evolution 
of a large number of particles the use of an efficient method for the evaluation of forces is 
necessary. For our simulations we used a tree-like algorithm of the kind firstly proposed by 
Barnes ct al. O]. In this algorithm, the force acting on a particle due to its neighbors is 
evaluated exactly, while that due to distant particles is evaluated by means of a multipolar 
expansion. 

The nature of the gravitational potential requires a time integrator, which is both flexible 
and efficient. Actually, the short range divergence can produce fast changes in the dynamical 
variables of a particle, whose motion should be therefore followed with a very short time step. 
Nevertheless, most of the particles may have a regular and smooth motion which doesn't need 
such a careful integration. The solution adopted is thus an algorithm, based on a variable 
and individual time-step. In order to avoid time steps approaching to zero when a very 
close encounter between two particles occurs, a smoothing in the potential is introduced at 
very small scales. This, often-used, technique introduces of course an artificial length scale 
in the simulations. In order to save the Newtonian behavior of the potential, such a scale 
must be taken much smaller than the average distance between nearest neighbors A. In the 
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simulations presented here, the smoothing scale is approximately 1/32 A. We have checked on 
some simulations with 2 12 particles that using a smaller smoothing length does not influence 
our results; however it slows down even more the code. 

The simulations presented in this paper were performed using the tree-code described 
in [H| which we adapted to represent an infinite system by means of periodic boundary 
conditions according to Ewald prescription (see e.g. flij ). 

A small smoothing scale and the absence of cosmological expansion, together with the 
requirement of high accuracy imply computations heavier than the usual cosmological ones. 
This is because we consider the generation of structures from a strongly non linear dynamics 
starting from random initial positions. This will limit the number of particles N in our systems, 
but all our results will be analyzed with respect to the dependence on N. Our simulations 
are presented using arbitrary units for mass (m u ), length (r u ) and time (t u ). The Newtonian 
equations depend on a single constant with physical dimensions (G), therefore we can set 
arbitrarily two of the units and the third is uniquely defined according to r^t^m^ 1 = G. 

The left column of fig. |l| shows the projection onto the xy-plane of typical snapshots of 
a system with N = 2 15 particles at different times. All our simulations are performed with 
a constant number density no = N/L 3 = lr~ 3 . With this choice, a typical crossing time 
r = (mnoG) -1 / 2 corresponds to lt u . On the right column we plot the conditional average 
density T*(r, t) of the system at the same time of the corresponding snapshot. T*(r, t) is 
the average density in spheres of radius r centered on a particle i of the system. It can be 
defined as: 

r * (r ' t] = w^n E E °( r - r *(*) + in i) ( 3 ) 

where O(r) = 47rr 3 /3 and r<(t), rj(t) are the positions at time t of particles i and j respectively. 

In the eq. |^, particles i and j belong to the simulation box, and the summation on the 
integer vectors n allow to take into account the replicas of the particle j. The three lines in 
fig. correspond to simulations with N = 2 9 (red), N = 2 12 (blue), N = 2 15 (black). Initial 
conditions (a) are given by N point masses located randomly and at rest into the simulation 
box. In the initial phase (b) the system develops clustered structures on small scales. Cor- 
respondingly the r*(r, t) acquires amplitude below the homogeneity scale ro, defined by the 
condition r*(?- , t) — 2n . The scale r can be seen as a typical clustering scale; the reason for 
such definition will be clearer with the definition of £(r, <).The red circle in the correspond- 
ing snapshot has radius ro, and roughly identifies the scale at which the fluctuations in the 
number density are of the same order of the average value uq. 

As long as the cluster sizes are much smaller than the size of the box they continue merging 
and forming bigger and bigger structures (c). Accordingly the homogeneity scale r increases. 
Eventually, the typical size of the few remaining aggregates becomes comparable with the box 
size and they collapse into a single big cluster which contains almost all the matter (d). In 
order to give a rough estimate the time scale for the formation of the final cluster, we can define 
a transition time tf, such that ro(i/) = L/4. For simulations with the same mass density, 
larger N corresponds to larger box size L; the final cluster is therefore larger and tf is greater. 
After the time tf the system reaches a genuine statistically-stationary state. The accuracy of 
our simulations can be checked by energy conservation. It is a bit tricky to evaluate energy 
conservation, because in an infinite system there is no absolute value for the potential. Then 
the relative error in energy AE/E is not a meaningful quantity. For this reason we use the 
quantity AE(t)/ AK(t), the ratio between the error in the total energy and the increase in 
kinetic energy. Typically, when the final cluster forms the energy conservation is within 1% . 

In the first three snapshots of fig. |l| the measured T*(r, t) function does not depend on 
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Fig. 1 - Left: projection onto the xy-p\ane at different times of a system with N = 2 particles in a 
simulation box with side L = 32. (a) t = 0; (b) t = 0.6; (c) t = 1; (d) t = 3.4. Only 1/4 of the points 
are shown for the sake of clarity. Right: Conditional average density of systems with N = 2 9 (red), 
N — 2 12 (blue), iV = 2 15 (black). Arrows mark the "homogeneity scale". Note that the structures of 
the first three figures are size independent, while the size of the system influences the final cluster. 



the number N of particles in the simulation box. This is because the characteristic size of 
the structures is much smaller than the box size. On the other hand, we can see in the last 
panel (d) that the T* (r, t) corresponding to systems with different number of particles do not 
overlap any more. Therefore, the size of the final cluster depends on the box size. 

In figure |^ we show the kinetic energy per particle K(t) versus time for the three simu- 
lations. We have marked the values of the kinetic energy corresponding to the snapshots in 
fig. [I]. Again, we can distinguish two regimes in the time evolution of K(t): (i) during the 
initial collapse phase K(t) show the same growth for the three simulations while (ii) at a later 
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Fig. 2 - Kinetic energy per particle for simulations with N = 2 9 (red), N = 2 12 (blue), N = 2 15 
(black). Circles identify the times used in fig. |l| Note that the dynamics of the simulations becomes 
different only at large times, when the cluster size is comparable to the box size. 



time, the equilibrium value depends strongly on the number N of particles. The breakdown 
of the Gibbs statistical treatment is reflected here by the fact that equilibrium properties of 
the system, like K(t — > oo) and T*(r, t — ► oo), do not posses a well-defined thermodynamic 
limit. In other words, they do not posses a finite limit for N — > oo and L — ► oo with N/L 3 
kept constant. On the other hand, the same quantities appear to be independent on the size 
of the system for times smaller than tf. This suggests that the transient phase, during which 
collapse occurs, does posses properties with a well behaved limit in the infinite system. This 
implies that the transient clustering is a well defined statistical properties of the system. 

In order to characterize the very small density fluctuations at the largest scales of the 
simulations, we make use the integrated two-point correlation function of density fluctuations 
£(r, t) — r*(r, t)/no — 1. It measures the fluctuations above the average value no inside a 
sphere of radius r. Our definition of rg then corresponds to |£(ro, t)\ = 1, which is a standard 



definition for the homogeneity scale 17 In fig. |3| we show the absolute value |£(r, t)\ versus r 
for the system with N = 2 15 . The lowest (black) line represents the correlation function of the 
initial configuration, which scales as r -3 / 2 as expected for points distributed randomly. The 
other lines corresponds to times t spaced by 1/15 t u . At the largest scales of our simulations, 
say for r > 10r u , the time evolution of the correlation function appears to be well described 
by the linear approximation for the equation of a self-gravitating fluid ]l6[ |, which predicts: 

f(r, t) = £(r, 0) * cosh 2 (t/r/) (4) 



with T[ ps (inGpo) ~? = l/y4wt u . Eq. ^ describes the linear regime in which the time 
evolution of the system consists solely in an exponential amplification of the initial density 
fluctuations. Such behavior is indicated in the figure, by the upward arrow labeled by F. 
It is worth noting that in this regime, overdensities increase while underdensities actually 
decrease; however, in fig. || we have plotted the absolute value of £(r, t), and therefore we see 
an enhancement of all the fluctuations. 

At small scales the dynamics of the system shows completely different properties. For an 
initial unclustered random distribution the discrete character of points of the system dominates 
the clustering process at small scales. At the very beginning, i.e. for t < i„/4, the fluctuations 
grow very fast due to the fall of each particle in the direction of the nearest particle. Actually, 
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Fig. 3 - Evolution of |£(r)| between t = and t = 1, spaced by time intervals At = 1/15 for a system 
with TV = 2 15 particles. Arrow P identifies the initial phase dominated by two particles interactions; 
F is associated with fluid-like evolution at large scales, but the main dynamics for cluster formation is 
due to the process indicated by G, which effectively shifts the granularity from small scales to large 
scales. 



Holtsmark distribution shows that the force acting on a particle in an infinite Poisson system 
is mainly due to its nearest neighbour. In this phase the dynamics is therefore dominated by 
two-body interactions. For times greater than approximately t u /A and at intermediate scales, 
the r*(r, t) function satisfies the following equation: 

r*(r, t+ At) =r*(e"^Y, t). (5) 

with a measured time-scale r g = 0.44 t u , and a similar relation also holds for £(r, t). This 
behavior corresponds to a rigid translation along the r axis in a logarithmic plot. This indicates 
a sort of similarity of structures of different size at different times. Although it is different 
from the self-similarity as usually intended, it is a clear indication of a sort of scale invariance 
in the dynamics. Eq. (|5|) states that the correlations developed by aggregates of size r at time 

-At 

t + At are the same as those developed by aggregates of size re T9 at time t, and that the 
time evolution of such correlations is the same. 

Visual inspection of the evolution of the system (fig. |l|) suggests that clusters formed at a 
given time are the new discrete elements for the iteration of the clustering process at a larger 
scale. We propose a simple argument which can account for the behavior described by Eq. (Q) . 
Consider a pair of particles at rest with mass m u and at distance r u , as in the initial condition 
of our simulations. Since the leading contribution of the force acting on a particle is mainly 
due to its nearest neighbor |l8fl, the time scale tq for their relative motion is: 




After At oc t most of the particles are much closer to their nearest neighbor. We will 
consider such pairs as "bound systems" approximately at rest and distributed randomly with 
typical distance n « 2 1 / 3 r u and total mass mj = 2m u . We can iterate the same scheme, 
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treating the pairs as the new "particles" of the system. The corresponding time scale T\ for 
the relative motion is: 

n = y/ri/Grm = v / 2r3/G2m„ = r . (7) 

Iterating the process we find a sequence of "bound systems" with increasing mass m n = 
2 n m u and length-scale r„ = 2™/ 3 r M , all with the same time-scale r„ = r . Therefore, being 
n ~ i/ T 0i the scale at which this activity takes place r(t) will obey the following equation: 

r(t) cx 2^r u (8) 

which accounts for the exponential dependence of Eq. (|J) . 

In summary, from the present study we can draw the following conclusions: (i) we address 
the issue of the cluster formation in a gravity driven dynamics in a system with random 
initial positions and periodic boundary conditions. Strictly speaking, the system has no 
well defined thermodynamic limit at asymptotic times. However, it shows a well defined, 
thermodynamically stable clustering in the transient regime; (ii) we identify the mechanisms 
of structure formation: initially the system tends to form pairs These pairs then interact with 
each other and gradually form larger and larger structures, shifting the granular structure 
from small to large scales. This process appears to be finally dominating over the fluid like 
dynamics, which tends to amplify the initial large scale fluctuations. These results point to 
a fundamental importance of the granular mechanism, which appears to lead to a dynamics 
with elements of self similarity. 

To check the validity of our results, we have also repeated some of the simulations with 
the GADGET code @. 

* * * 

This work was partially supported by the INFM under the project Clustering and by the 
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